options nofmterr; 
libname meps05 'E:\MEPS\new_basevar_2005'; 
libname para 'E:\MEPS\Final\para';

proc format;
    value female   1 = "a Women" 2= "b men";
    value race     1 = "a White" 2="b other";
    value metro    1 = "a metro" 2 = "b nonmetro";
    value marital  1 = "a Married" 2 = "b other";
    value yesno    1 = "a Yes" 2="b No";
    value occ      1 = "a metro" 2 = "b non-metro ";
	value region   1 = "a northeast" 2 = "b midwest" 3 = "c south" 4 = "d west" 5 = "e missing"; 
	value famsize  1 = "a 1-2" 2 = "b 3-4" 3 = "c 5+"; 
    value age      1 = "a < 40 yrs" 2 = "b 40 - 49" 3 = "c 50 - 64" ;
    value chronic  0 = "a no chronic disease" 1 = "b 1 chronic disease" 2 = "c 2 or more";
	value income   0 = "a <20,000" 1 = "b 20,000-30,000" 2 = "c 30,000-50,000"  3 ="d 50,000+";
    value educ     1 = "a <HS" 2 = "b HS" 3 = "c some college + ";
    value ind      1 = "a <75%" 2 = "b 75%-90%" 3 = "90% + ";


data obese; 
set meps05.new_basevar_2005; 
obese_pop = 0; 
if (bmi =4) then obese_pop = 1; 
else obese_pop = 0; 
if (bp = 1 ! heart = 1 ! diabetes = 1 ! osteo_arth = 1 ! 
asthma = 1 ) then chronic_pop=1; 
else chronic_pop=0; 
length dep_anx_sz 3.; label dep_anx_sz = "Any dep,anx,sz"; 
if (dep = 1 ! anxiety = 1 ! schiz = 1) then dep_anx_sz = 1; 
else dep_anx_sz = 0; 
if obese_pop = 1 & (inscov ne 2) & INSCOV>=0 & eligpop =1 & PERWT>0 & 21<age<65 
& tot_exp>=0 & female>=0 & race>=0 & marr>=0 & educyr>=0 & metro>=0 & income>=0 
& CURR_SMK>=0 & exercise>=0 & dep_anx_sz>=0 & alzhmr>=0 & asthma>=0 & arthritis>=0 & cancer>=0 & 
emphysema>=0 & diabetes>=0 & heart>=0 & bp>=0 & osteo>=0 & stroke>=0 
& employ>=1 & occcat>=0 & indcat>=0; 
if inscov=3 then inscov=0; 
visit=(tot_exp>0); 
female=2-female; 
metro=2-metro; 
exercise=2-exercise; 
CURR_SMK=2-CURR_SMK; 
employ=2-employ; 
*************************** 
* # of comorbid conditions* 
***************************; 
comorb_nbr = sum (alzhmr,asthma,arthritis,cancer,emphysema,diabetes,heart,bp,osteo,stroke); 
if (comorb_nbr = .) then comorb_nbr = 0; 
run; 
data try; 
set obese; 
if indcat not in (14,15); 
if indcat=1 then ind_cnts=56/(56+35); 
if indcat in (2,4) then ind_cnts=82/(82+14); 
if indcat=3 then ind_cnts=59/(59+37); 
if indcat=5 then ind_cnts=72/(72+20); 
if indcat=6 then ind_cnts=76/(76+19); 
if indcat in (7,10) then ind_cnts=88/(88+8); 
if indcat=8 then ind_cnts=86/(86+10); 
if indcat=9 then ind_cnts=74/(74+20); 
if indcat in (11,12) then ind_cnts=59/(59+32); 
if indcat=13 then ind_cnts=91/(91+5); 
ind_cnts=ind_cnts; 
if occcat in (1,2) then occ=1; /* management, professional */ 
else occ=0; 
if marr=1 then marital=1; /* married */ 
else marital=0; 
if race=1 then white=1; else white=0; 
if race=2 then black=1; else black=0; 
if race=3 then hispc=1; else hispc=0; 
if race=4 then otherace=1; else otherace=0; 
ln_exp=log(tot_exp+1); 
tot_exp=tot_exp/1000; 
ln_income= log(income+1); 
income=income/1000; 
age2=age*age/100; 
if metro=1 then metreg=region; 
else metreg=4+region; 
if region=1 then reg_ne=1; else reg_ne=0; 
if region=2 then reg_mw=1; else reg_mw=0; 
if region=3 then reg_so=1; else reg_so=0; 
run; 
data add; 
set meps05.h97 (keep= DUPERSID FAMSZEYR); 
run; 
data try; 
merge add try(in=A); 
by DUPERSID;
if A; 
run; 
proc freq data=try; 
tables FAMSZEYR; 
run; 

data meps05.final;
set try;
keep age age2 comorb_nbr dep_anx_sz female white ln_income CURR_SMK educyr marital FAMSZEYR 
	reg_ne reg_mw reg_so ind_cnts occ inscov visit ln_exp; 
run;

/*proc corr data=meps05.final;
var marital FAMSZEYR;
run;

data temp;
set meps05.final;
if FAMSZEYR>1 then married=1; else married=0;
proc freq data=temp;
tables marital*married;
run;
proc corr data=temp;
var marital married;
run;
run;*/


data desc; 
set try; 
********************* 
* Age grouped * 
*********************; 
if (0 <= age <= 39) then age_grp = 1; 
else if (40 <= age <= 49) then age_grp = 2; 
else if (50 <= age <= 64) then age_grp = 3; 
************************* 
* education * 
*************************; 
if ( 0 <= educyr <= 11) then educ_grp = 1; /* LT HS */ 
else if ( educyr = 12) then educ_grp = 2; /* HS */ 
else if ( 13 <= educyr <= 17) then educ_grp = 3; /* some College + */ 
************************* 
* income * 
*************************; 
if ( 0 <= income < 20) then income_grp = 0; /* LT 20,000 */ 
else if ( 20 <= income < 30) then income_grp = 1; /* LT 30,000 */ 
else if ( 30 <= income < 50) then income_grp = 2; /* LT 50,000 */ 
else if ( 50 <= income ) then income_grp = 3; /* 50,000 + */ 
************************* 
* comorbidities * 
*************************; 
if ( comorb_nbr = 0) then comorb_grp = 0; /* No comorbidities */ 
else if ( comorb_nbr = 1) then comorb_grp = 1; /* 1 */ 
else if ( comorb_nbr >= 2) then comorb_grp = 2; /* 2+ */ 
************************* 
* industry insurance * 
*************************; 
if ( ind_cnts < 0.75) then ind_grp = 0; /* <75% */ 
else if ( 0.75 <= ind_cnts < 0.9) then ind_grp = 1; /* 75%-90% */ 
else if ( ind_cnts >= 0.9) then ind_grp = 2; /* 90%+ */ 
************************* 
* family size * 
*************************; 
if ( 1 <= FAMSZEYR <= 2) then fam_grp = 1; /* 1-2 */ 
else if ( 3 <= FAMSZEYR <= 4) then fam_grp = 2; /* 3-4 */ 
else if ( FAMSZEYR >= 5) then fam_grp = 3; /* 5+ */ 
********************* 
* Recode 0s to 2 * 
*******************; 
array reco inscov visit dep_anx_sz exercise curr_smk female white metro marital occ; 
do over reco; 
if (reco = 0) then reco = 2; 
end; 
length tot 3.; label tot = "1 for everyone"; 
tot = 1; 
run; 
************************** 
* Base dataset to append * 
**************************; 
data all; 
count = .; percent =.; 
format count f9.0 percent 5.1; 
%macro frq(var); 
proc freq data=desc; 
tables &var /out=d ; 
format inscov visit dep_anx_sz exercise curr_smk yesno. metro metro. region region. occ occ. 
female female. white race. marital marital. comorb_grp chronic. 
age_grp age. educ_grp educ. income_grp income. ind_grp ind.; 
title "Describe study population obese observations"; 
*********************** 
* APpend dataset * 
***********************; 
data d; set d; keep count percent; 
proc append base = all data = d; 
%mend; 
%frq(tot); 
%frq(inscov); 
%frq(visit); 
%frq(educ_grp); 
%frq(age_grp); 
%frq(income_grp); 
%frq(female); 
%frq(white); 
%frq(comorb_grp); 
%frq(dep_anx_sz); 
%frq(exercise); 
%frq(CURR_SMK); 
%frq(marital); 
%frq(metro); 
%frq(fam_grp); 
%frq(region); 
%frq(metreg); 
%frq(ind_grp); 
%frq(occ); 
****************** 
* Export results * 
******************; 
proc export data=all 
outfile='E:\MEPS\new_basevar_2005\desc.dbf' replace dbms=dbf; 
run; 
************************** 
* Base dataset to append * 
**************************; 
data all; 
count1 = .; rowper1 =.; 
count2 = .; rowper2 = .; 
chisq = .; pvalue = .; 
length sig $3.; 
totn = .; totper = .; 
format count1 count2 totn 6.0 
rowper1 rowper2 totper 8.1; 
run; 
%macro crstab(var); 
proc freq data=desc; 
tables &var * inscov /chisq out=d outpct; 
output out=stat chisq; 
format dep_anx_sz exercise curr_smk yesno. metro metro. region region. occ occ. 
female female. white race. marital marital. comorb_grp chronic. 
age_grp age. educ_grp educ. income_grp income. ind_grp ind.; 
title "Crosstab with Insurance"; 
*************************** 
* Rename variables for * 
* Easy tabulation * 
***************************; 
data one (rename =(count=count1 pct_row=rowper1)) 
two (rename =(count=count2 pct_row=rowper2)); 
set d (keep = count pct_row inscov); 
if (inscov = 1) then output one; 
if (inscov = 2) then output two; 
********************* 
* Get total % * 
*********************; 
data total; set d (keep = count percent &var); by &var; 
if (first.&var) then do; 
totn = 0; totper = 0; 
end; 
retain totn totper; 
totn = count + totn; 
totper = Percent + totper; 
if (last.&var) then do; 
keep totn totper; 
output; 
end; 
*********************** 
* Only desired stat * 
***********************; 
data stat; 
set stat (keep = _pchi_ p_pchi rename = (p_pchi = pvalue _pchi_ = chisq)); 
length sig $3.; 
if (0 <= pvalue < .001) then sig = "***"; 
else if (.001 <= pvalue < .01) then sig = "**"; 
else if (.01 <= pvalue < .05) then sig = "*"; 
data total; merge one two total; 
data total; set stat total; 
format count1 count2 totn 6.0 
rowper1 rowper2 totper 8.1; 
keep count1 rowper1 count2 rowper2 
totn totper chisq pvalue sig; 
proc append base=all data=total; 
%mend crstab; 
%crstab(inscov); 
%crstab(visit); 
%crstab(educ_grp); 
%crstab(age_grp); 
%crstab(income_grp); 
%crstab(female); 
%crstab(white); 
%crstab(comorb_grp); 
%crstab(dep_anx_sz); 
%crstab(exercise); 
%crstab(CURR_SMK); 
%crstab(marital); 
%crstab(fam_grp); 
%crstab(metro); 
%crstab(region); 
%crstab(metreg); 
%crstab(ind_grp); 
%crstab(occ); 
****************** 
* Export results * 
******************; 
proc export data=all 
outfile='E:\MEPS\new_basevar_2005\crossinscov.dbf' replace dbms=dbf; 
run; 
************************** 
* Base dataset to append * 
**************************; 
data all; 
count1 = .; rowper1 =.; 
count2 = .; rowper2 = .; 
chisq = .; pvalue = .; 
length sig $3.; 
totn = .; totper = .; 
format count1 count2 totn 6.0 
rowper1 rowper2 totper 8.1; 
run; 
%macro crstab(var); 
proc freq data=desc; 
tables &var * visit /chisq out=d outpct; 
output out=stat chisq; 
format dep_anx_sz exercise curr_smk yesno. metro metro. occ occ. 
female female. white race. marital marital. comorb_grp chronic. 
age_grp age. educ_grp educ. income_grp income. ind_grp ind.; 
title "Crosstab with Visit"; 
*************************** 
* Rename variables for * 
* Easy tabulation * 
***************************; 
data one (rename =(count=count1 pct_row=rowper1)) 
two (rename =(count=count2 pct_row=rowper2)); 
set d (keep = count pct_row visit); 
if (visit = 1) then output one; 
if (visit = 2) then output two; 
********************* 
* Get total % * 
*********************; 
data total; set d (keep = count percent &var); by &var; 
if (first.&var) then do; 
totn = 0; totper = 0; 
end; 
retain totn totper; 
totn = count + totn; 
totper = Percent + totper; 
if (last.&var) then do; 
keep totn totper; 
output; 
end; 
*********************** 
* Only desired stat * 
***********************; 
data stat; 
set stat (keep = _pchi_ p_pchi rename = (p_pchi = pvalue _pchi_ = chisq)); 
length sig $3.; 
if (0 <= pvalue < .001) then sig = "***"; 
else if (.001 <= pvalue < .01) then sig = "**"; 
else if (.01 <= pvalue < .05) then sig = "*"; 
data total; merge one two total; 
data total; set stat total; 
format count1 count2 totn 6.0 
rowper1 rowper2 totper 8.1; 
keep count1 rowper1 count2 rowper2 
totn totper chisq pvalue sig; 
proc append base=all data=total; 
%mend crstab; 
%crstab(inscov); 
%crstab(visit); 
%crstab(educ_grp); 
%crstab(age_grp); 
%crstab(income_grp); 
%crstab(female); 
%crstab(white); 
%crstab(comorb_grp); 
%crstab(dep_anx_sz); 
%crstab(exercise); 
%crstab(CURR_SMK); 
%crstab(marital); 
%crstab(fam_grp); 
%crstab(metro); 
%crstab(region); 
%crstab(metreg); 
%crstab(ind_grp); 
%crstab(occ); 
****************** 
* Export results * 
******************; 
proc export data=all 
outfile='E:\MEPS\new_basevar_2005\crossvisit.dbf' replace dbms=dbf; 
run; 


proc qlim data=try OUTEST=bivest COVOUT; 
model inscov = age age2 comorb_nbr dep_anx_sz female white ln_income CURR_SMK educyr marital FAMSZEYR reg_ne reg_mw reg_so ind_cnts occ; 
model visit = age age2 comorb_nbr dep_anx_sz female white ln_income CURR_SMK educyr marital FAMSZEYR reg_ne reg_mw reg_so inscov; 
endogenous inscov visit ~ discrete; 
	 INIT  _Rho=-0.05831;
class female white metro marital occ region reg_ne reg_mw reg_so 
CURR_SMK exercise dep_anx_sz; 
output out=resultlnln marginal XBETA proball; 
run;


data parm;
set bivest;
if _type_='PARM';
run;

data testrun;
  if _N_ = 1 then set parm(drop = _NAME_ _TYPE_ _STATUS_);
set resultlnln;
run;

data meps05.para_margin; 
set testrun;

/* inscov equation */
XB_age_2=XBETA_inscov+   inscov_age - age2*(inscov_age2) + (age+1)*(age+1)*(inscov_age2)/100;
mar_inscov_age=cdf('normal',XB_age_2)-cdf('normal',XBETA_inscov);

XB_comorb_1=XBETA_inscov+  inscov_comorb_nbr;
mar_inscov_comorb=cdf('normal',XB_comorb_1)-cdf('normal',XBETA_inscov);

XB_dep_0=XBETA_inscov-dep_anx_sz*( inscov_dep_anx_sz);
XB_dep_1=XBETA_inscov+(1-dep_anx_sz)*( inscov_dep_anx_sz);
mar_inscov_dep=cdf('normal',XB_dep_1)-cdf('normal',XB_dep_0);

XB_female_0=XBETA_inscov-female*(inscov_female);
XB_female_1=XBETA_inscov+(1-female)*(inscov_female);
mar_inscov_female=cdf('normal',XB_female_1)-cdf('normal',XB_female_0);

XB_white_0=XBETA_inscov-white*( inscov_white);
XB_white_1=XBETA_inscov+(1-white)*( inscov_white);
mar_inscov_white=cdf('normal',XB_white_1)-cdf('normal',XB_white_0);

XB_ln_income_1=XBETA_inscov+ 0.1*  inscov_ln_income; /* increase income by 10% */
mar_inscov_ln_income=cdf('normal',XB_ln_income_1)-cdf('normal',XBETA_inscov);

XB_curr_smk_0=XBETA_inscov-curr_smk*( inscov_curr_smk);
XB_curr_smk_1=XBETA_inscov+(1-curr_smk)*(inscov_curr_smk);
mar_inscov_curr_smk=cdf('normal',XB_curr_smk_1)-cdf('normal',XB_curr_smk_0);

XB_educ_1=XBETA_inscov+  inscov_educyr;
mar_inscov_educ=cdf('normal',XB_educ_1)-cdf('normal',XBETA_inscov);

XB_marital_0=XBETA_inscov-marital*( inscov_marital);
XB_marital_1=XBETA_inscov+(1-marital)*( inscov_marital);
mar_inscov_marital=cdf('normal',XB_marital_1)-cdf('normal',XB_marital_0);

XB_FAMSZEYR_1=XBETA_inscov+  inscov_FAMSZEYR;
mar_inscov_FAMSZEYR=cdf('normal',XB_FAMSZEYR_1)-cdf('normal',XBETA_inscov);

XB_reg_ne_0=XBETA_inscov-reg_ne*( inscov_reg_ne);
XB_reg_ne_1=XBETA_inscov+(1-reg_ne)*(inscov_reg_ne);
mar_inscov_reg_ne=cdf('normal',XB_reg_ne_1)-cdf('normal',XB_reg_ne_0);

XB_reg_mw_0=XBETA_inscov-reg_mw*( inscov_reg_mw);
XB_reg_mw_1=XBETA_inscov+(1-reg_mw)*(inscov_reg_mw);
mar_inscov_reg_mw=cdf('normal',XB_reg_mw_1)-cdf('normal',XB_reg_mw_0);

XB_reg_so_0=XBETA_inscov-reg_so*( inscov_reg_so);
XB_reg_so_1=XBETA_inscov+(1-reg_so)*(inscov_reg_so);
mar_inscov_reg_so=cdf('normal',XB_reg_so_1)-cdf('normal',XB_reg_so_0);

XB_ind_1=XBETA_inscov+  inscov_ind_cnts * 0.05; /* assuming increasing 5% in industry coverage */
mar_inscov_ind=cdf('normal',XB_ind_1)-cdf('normal',XBETA_inscov);

XB_occ_0=XBETA_inscov-occ*( inscov_occ);
XB_occ_1=XBETA_inscov+(1-occ)*( inscov_occ);
mar_inscov_occ=cdf('normal',XB_occ_1)-cdf('normal',XB_occ_0);

/* visit equation */
XB_age_1=XBETA_visit + visit_age - age2*(visit_age2) + (age+1)*(age+1)*(visit_age2)/100;
mar_visit_age=cdf('normal',XB_age_1)-cdf('normal',XBETA_visit);

XB_comorb_1=XBETA_visit+  visit_comorb_nbr;
mar_visit_comorb=cdf('normal',XB_comorb_1)-cdf('normal',XBETA_visit);

XB_dep_0=XBETA_visit-dep_anx_sz*( visit_dep_anx_sz);
XB_dep_1=XBETA_visit+(1-dep_anx_sz)*( visit_dep_anx_sz);
mar_visit_dep=cdf('normal',XB_dep_1)-cdf('normal',XB_dep_0);

XB_female_0=XBETA_visit-female*( visit_female);
XB_female_1=XBETA_visit+(1-female)*(visit_female);
mar_visit_female=cdf('normal',XB_female_1)-cdf('normal',XB_female_0);

XB_white_0=XBETA_visit-white*(  visit_white);
XB_white_1=XBETA_visit+(1-white)*(  visit_white);
mar_visit_white=cdf('normal',XB_white_1)-cdf('normal',XB_white_0);

XB_ln_income_1=XBETA_visit + 0.1* visit_ln_income; /* increase income by 10% */
mar_visit_ln_income=cdf('normal',XB_ln_income_1)-cdf('normal',XBETA_visit);

XB_curr_smk_0=XBETA_visit-curr_smk*( visit_curr_smk);
XB_curr_smk_1=XBETA_visit+(1-curr_smk)*(visit_curr_smk);
mar_visit_curr_smk=cdf('normal',XB_curr_smk_1)-cdf('normal',XB_curr_smk_0);

XB_educ_1=XBETA_visit+ visit_educyr;
mar_visit_educ=cdf('normal',XB_educ_1)-cdf('normal',XBETA_visit);

XB_marital_0=XBETA_visit-marital*( visit_marital);
XB_marital_1=XBETA_visit+(1-marital)*( visit_marital);
mar_visit_marital=cdf('normal',XB_marital_1)-cdf('normal',XB_marital_0);

XB_FAMSZEYR_1=XBETA_visit+  visit_FAMSZEYR;
mar_visit_FAMSZEYR=cdf('normal',XB_FAMSZEYR_1)-cdf('normal',XBETA_visit);

XB_reg_ne_0=XBETA_visit-reg_ne*( visit_reg_ne);
XB_reg_ne_1=XBETA_visit+(1-reg_ne)*(visit_reg_ne);
mar_visit_reg_ne=cdf('normal',XB_reg_ne_1)-cdf('normal',XB_reg_ne_0);

XB_reg_mw_0=XBETA_visit-reg_mw*( visit_reg_mw);
XB_reg_mw_1=XBETA_visit+(1-reg_mw)*(visit_reg_mw);
mar_visit_reg_mw=cdf('normal',XB_reg_mw_1)-cdf('normal',XB_reg_mw_0);

XB_reg_so_0=XBETA_visit-reg_so*( visit_reg_so);
XB_reg_so_1=XBETA_visit+(1-reg_so)*(visit_reg_so);
mar_visit_reg_so=cdf('normal',XB_reg_so_1)-cdf('normal',XB_reg_so_0);

XB_inscov_0=XBETA_visit-inscov*( visit_inscov);
XB_inscov_1=XBETA_visit+(1-inscov)*( visit_inscov);
mar_visit_inscov=cdf('normal',XB_inscov_1)-cdf('normal',XB_inscov_0);

run;


proc means data=meps05.para_margin; 
var mar_inscov_age mar_inscov_comorb mar_inscov_dep mar_inscov_female mar_inscov_white mar_inscov_ln_income mar_inscov_curr_smk
    mar_inscov_educ mar_inscov_marital mar_inscov_FAMSZEYR mar_inscov_reg_ne mar_inscov_reg_mw mar_inscov_reg_so mar_inscov_ind mar_inscov_occ
	mar_visit_age mar_visit_comorb mar_visit_dep mar_visit_female mar_visit_white mar_visit_ln_income mar_visit_curr_smk
    mar_visit_educ mar_visit_marital mar_visit_FAMSZEYR mar_visit_reg_ne mar_visit_reg_mw mar_visit_reg_so mar_visit_inscov;
run;

data expenditure;                                                                                                                    
set testrun;                                                                                                                     
if visit>0;                                                                                                                             
Prob2_inscov=cdf('normal',XBETA_inscov);                                                                                                
Prob1_inscov=1-Prob2_inscov;                                                                                                            
A1=XBETA_visit-visit_inscov*inscov+visit_inscov; /* XA*bA+thetaA*/                                                                                
A0=XBETA_visit-visit_inscov*inscov; /* XA*bA*/                                                                                               
Z1=pdf('normal',A1)/cdf('normal',A1);                                                                                                   
Z2=pdf('normal',XBETA_inscov)/cdf('normal',XBETA_inscov);                                                                               
Z3=pdf('normal',A0)/cdf('normal',A0);                                                                                                   
Z4=-pdf('normal',XBETA_inscov)/(1-cdf('normal',XBETA_inscov));                                                                          
S=Prob2_inscov*pdf('normal',A1)/cdf('normal',A1) +                                                                                      
Prob1_inscov*pdf('normal',A0)/cdf('normal',A0);                                                                                         
run;                                                                                                                                    
                                                                                                                                        
                                                                                                                                        
data ex1 ex0;                                                                                                                           
set expenditure;                                                                                                                        
if inscov=1 then output ex1;                                                                                                            
else output ex0;                                                                                                                        
run;                                                                                                                                    
                                                                                                                                        
data ex1;                                                                                                                               
set ex1(rename= (Z1=ZZ1 Z2=ZZ2));                                                                                                       
                                                                                                                                        
data ex0;                                                                                                                               
set ex0(rename= (Z3=ZZ1 Z4=ZZ2));                                                                                                       
run;                                                                                                                                    
                                                                                                                                        
data combine;                                                                                                                           
set ex1(keep= ln_exp inscov age age2 comorb_nbr dep_anx_sz female                                                                       
white ln_income curr_smk educyr marital FAMSZEYR ind_cnts occ 
reg_ne reg_mw reg_so ZZ1 ZZ2)                                                                                                                              
	ex0(keep= ln_exp inscov age age2 comorb_nbr dep_anx_sz female                                                                   
white ln_income curr_smk educyr marital FAMSZEYR ind_cnts occ 
reg_ne reg_mw reg_so ZZ1 ZZ2);                                                                   
run;                                                                                                                                    
                                                                                                                                        
proc corr data=combine;                                                                                                                 
var ZZ1 ZZ2;                                                                                                                            
run;                                                                                                                                    
                                                                                                                                        
proc glm        data=combine;                                                                                                           
 title 'Proc glm Analysis';                                                                                                             
 model ln_exp = age age2 comorb_nbr dep_anx_sz female white ln_income                                                                   
CURR_SMK educyr inscov ZZ1 ZZ2;                                                                                                       
 output       out=fit p=yhat r=resid DFFITS=dfits COOKD=cookd;                                                                          
run;                                                                                                                                    
quit;

proc contents data=meps05.final;
run;

data meps05.testrun;
set testrun;
run;

data meps05.bivest;
set testrun;
run;


proc genmod    data=combine;                                                                                                           
 title 'Proc genmod Analysis';                                                                                                             
 model ln_exp = age age2 comorb_nbr dep_anx_sz female white ln_income                                                                   
CURR_SMK educyr inscov ZZ1 ZZ2 /dist = normal  link = identity CovB;                                                                                                       
run;                                                                                                                                    
quit;

proc means data=expenditure;
var age;
run;

proc gplot data=para.para_margin;
plot mar_inscov_age*age mar_inscov_comorb*comorb_nbr mar_inscov_ln_income*ln_income mar_inscov_educ*educyr
 mar_inscov_FAMSZEYR*FAMSZEYR mar_inscov_ind*ind_cnts;
run;
quit;

proc gplot data=para.para_margin;
plot mar_visit_age*age mar_visit_comorb*comorb_nbr mar_visit_ln_income*ln_income mar_visit_educ*educyr
 mar_visit_FAMSZEYR*FAMSZEYR;
run;
quit;

goptions reset=all;
title "Marginal Effect of Age on Insurance";
axis1 label=("age");
axis2 label=("marginal effect");

proc gplot data=para.para_margin;
plot mar_inscov_age*age/haxis=axis1 vaxis=axis2;
run;
quit;

goptions reset=all;
title "Marginal Effect of Age on Utilization";
axis1 label=("age");
axis2 label=("marginal effect");

proc gplot data=para.para_margin;
plot mar_visit_age*age/haxis=axis1 vaxis=axis2;
run;
quit;
